clear; close
d = 1.2;    % m
hb = 0.3;   % m
g = 9.86;   % m/s^2

syms ha v0 theta;

% Numerical caculaition & plot
theta = 0.1:(1/20*pi):(pi/2-0.1);
ha = 0.2;

v0 = sqrt( (g.*d.^2.*(1+tan(theta).^2) ) ./ ( 2.*(ha-hb + d*tan(theta)) ));
x = 0:0.1:1.5;
grid on
for i = 1:length(theta)
    rad = theta(i);
    v = v0(i);
    
%     hold on;
%     subplot(1,2,1);
%     y = ha + tan(rad).*x - 1/2.*g.*x.^2./(v.*cos(rad)).^2;
%     plot(x,y);
%     legend;

    hold on;
%     subplot(1,2,2); 
    y = ha + tan(rad).*x - 1/2.*g.*x.^2./v.^2.*(1+tan(rad).^2);
    plot(x,y);
%     legend;
end

cat(1, v0, theta.*180/pi)